# replication code for results in "Politicization and Repsonsiveness in Executive Agencies"
require(ggplot2)
require(stargazer)
require(interplot)
require(data.table)
require(multiwayvcov)
require(grid)
require(gtable)

# setwd("/Users/klowande/")

###### Main Text ######

# Table 2: Modeling Agency Responsiveness to Members of Congress
load('d1.Rds')
load('d2.Rds')
load('d3.Rds')
load('d4.Rds')
load('d5.Rds')
m1 <- lm(rt_log ~  presparty * polit_sesrq + selin_dim1 + selin_dim2 + rep + senate 
         + enacted_in_billions + staff_in_thousands + wl.est + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=d3)
m2 <- lm(rt_log ~ agen  + presparty * polit_sesrq + rep + senate 
         + enacted_in_billions + staff_in_thousands + wl.est + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=d2)
m3 <- lm(rt_log ~ agen + mc.id + presparty * polit_sesrq 
         + enacted_in_billions + staff_in_thousands + wl.est + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=d1)
m4 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
         + enacted_in_billions + staff_in_thousands + wl.est + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=d1)
m5 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
         + enacted_in_billions + staff_in_thousands + wl.est  + majchair_r2 + minchair_r2 
         + seniority + maj, data=d4)
m6 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
         + enacted_in_billions + staff_in_thousands + wl.est + majchair_r2 + minchair_r2 
         + seniority + maj, data=d5)

# output
# stargazer(m1,m2,m3,m4,m5,m6,digits=2,out.header=F,out='OUTPUT-FILE PATH')

# Figure 1: Conditional Impact of President's Party on Response Time

figvals <- interplot(m4,"presparty","polit_sesrq",plot=F,hist=T,xmin=0,xmax=1)
kd <- density(d1$polit_sesrq)
dens <- data.frame(kd$x,kd$y)
a <- ggplot(figvals, aes(polit_sesrq, coef)) + 
    geom_line(data=figvals) +
    geom_ribbon(data=figvals,aes(ymin=lb,ymax=ub),alpha=.3, color= 'black') +
    ylab("Marginal Effect of President's Party") + geom_hline(yintercept=0,linetype="dashed") +
    xlab('') + theme(axis.title=element_text(size=17),axis.text=element_text(size=13))
b <- ggplot(data=d1[which(d1$polit_sesrq<=1),], aes(x=polit_sesrq)) +
    geom_histogram(data=subset(d1,presparty == 1),colour='black',fill='#a51417',bins=12,alpha=.75) +
    geom_histogram(data=subset(d1,presparty == 0),colour='black',fill='white',bins=12,alpha=.5) +
    ylab('Moderator Frequency') + xlab('Politicization') + 
    scale_x_continuous(breaks=c(0,0.25,0.5,0.75,1)) + coord_cartesian(xlim = c(-.02,1.02)) +
    annotate("text", label = "Red = President's Party", x = 0.82, y = 1900, size = 5, colour = "black") +
    annotate("text", label = "White = Opposition", x = 0.775, y = 1600, size = 5, colour = "black") +
    theme(axis.title=element_text(size=17),axis.text=element_text(size=13))
# output
# pdf("OUTPUT FILE PATH",height=10,width=6)
# grid.newpage()
# grid.draw(rbind(ggplotGrob(a), ggplotGrob(b), size = 'last'))
# dev.off()

###### Appendix ######

# TABLE A1: Modeling Legislator Contact
load('dt.Rds')
m1 <- lm(d.policy ~ agen + contact_cong +
           senate + republican + maj + presparty + seniority + chair + ranking + member,
         data=dt)
m2 <- lm(d.casework ~ agen + contact_cong +
           senate + republican + maj + presparty + seniority + chair + ranking + member,
         data=dt)
m3 <- lm(d.total ~ agen + contact_cong +
           senate + republican + maj + presparty + seniority + chair + ranking + member,
         data=dt)
m4 <- lm(d.policy ~ agen + contact_cong + mc.id +
           maj + presparty + seniority + chair + ranking + member,
         data=dt)
m5 <- lm(d.casework ~ agen + contact_cong + mc.id +
           maj + presparty + seniority + chair + ranking + member,
         data=dt)
m6 <- lm(d.total ~ agen + contact_cong + mc.id +
           maj + presparty + seniority + chair + ranking + member,
         data=dt)
# output
#stargazer(m1,m2,m3,m4,m5,m6,digits=2,out.header=F,out='OUTPUT FILE PATH')

# FIGURE E2: Workload and Estimated Workload
load('agency-wl-fixed.Rds')
f <-ggplot(data=wl[which(wl$agen==4 & wl$day<=16000),]) +
  geom_line(aes(x=day, y=workload)) + geom_line(color='red',aes(x=day, y=wl.est)) +
  xlab('Day') + ylab('Workload')
# output
#pdf("OUTPUT FILE PATH",height=3,width=8)
#f
#dev.off()

# FIGURE E3: Agency Workload Over Time
load('agency_lwl.Rds')
ef2a <- ggplot(data=lwl, aes(x=day, y=workload, group=agen, colour=agen)) + geom_line() + facet_grid(agen ~.) +
  geom_vline(xintercept = 14260, linetype="dotted") + geom_vline(xintercept = 14977, linetype="dotted") +
  ylab("Workload") + xlab("Day") + scale_x_continuous(breaks = c(14260,14977), labels = c("(a)","(b)")) +
  theme(legend.position="none")
load('agency_hwl.Rds')
ef2b <- ggplot(data=hwl, aes(x=day, y=workload, group=agen, colour=agen)) + geom_line() + facet_grid(agen ~.) +
  geom_vline(xintercept = 14260, linetype="dotted") + geom_vline(xintercept = 14977, linetype="dotted") +
  ylab("Workload") + xlab("Day") + scale_x_continuous(breaks = c(14260,14977), labels = c("(a)","(b)")) +
  theme(legend.position="none")
# output
#pdf("OUTPUT FILE PATH",height=8,width=4)
#ef2a
#dev.off()
#pdf("OUTPUT FILE PATH",height=8,width=4)
#ef2b
#dev.off()

# TABLE F1: Agency/Legislator Controls Omitted from Table 2 (see above)

# FIGURE F1: Conditional Responsiveness Time by Agency (see Table 2 above)

# FIGURE F2: Conditional Responsiveness Time by Year (see Table 2 above)

# TABLE F2: Primary model w/Dept. of Energy Included
load('de.Rds')
h1 <- lm(rt_log ~ agen + mc.id + presparty * polit_sesrq 
         + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=de)
h2 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
         + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=de)
#stargazer(h1,h2,digits=2,out.header=F,out="OUTPUT FILE PATH")

# TABLE F3: Primary model w/logged Number of Appointees as alt. IV
load('dh.Rds')
h3 <- lm(rt_log ~ agen + mc.id + presparty * lntot
         + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=dh)
h4 <- lm(rt_log ~ agen + mc.id + yr + presparty * lntot
         + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
         + seniority + maj, data=dh)
#stargazer(h3,h4,digits=2,out.header=F,out="OUTPUT FILE PATH")

# TABLE F4: Primary model as Negbin
h5 <- glm.nb(responsetime ~ agen + mc.id + presparty * polit_sesrq 
             + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
             + seniority + maj, data=d1)
h6 <- glm.nb(responsetime ~ agen + mc.id + yr + presparty * polit_sesrq 
             + enacted_in_millions + staffq + workload + casework + majchair_r2 + minchair_r2 
             + seniority + maj, data=d1)
#stargazer(h5,h6,digits=2,out.header=F,out="OUTPUT FILE PATH")

# TABLE F5: Modeling Agency Responsiveness to Members of Congress (Restricted Sample)
load('dr.Rds')
load('dr1.Rds')
load('dr2.Rds')
mr1 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
          + enacted_in_billions + staff_in_thousands + wl.est + casework + majchair_r2 + minchair_r2 
          + seniority + maj, data=dr)
mr2 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
          + enacted_in_billions + staff_in_thousands + wl.est  + majchair_r2 + minchair_r2 
          + seniority + maj, data=dr1)
mr3 <- lm(rt_log ~ agen + mc.id + yr + presparty * polit_sesrq 
          + enacted_in_billions + staff_in_thousands + wl.est + majchair_r2 + minchair_r2 
          + seniority + maj, data=dr2)
# output
# stargazer(mr1,mr2,mr3,digits=2,out.header=F,out="OUTPUT FILE PATH")

# FIGURE G1: Raw Plot
# see Hainmueller, Mummolo, Xu (2017): https://github.com/xuyiqing/interflex
source('hmx.R')
vars <- c("agen","idno","yr","enacted_in_millions","staffq","workload","casework",
          "majchair_r2","minchair_r2","seniority","maj")
ff1 <- inter.rawplot(Y = "rt_log", D = "presparty", X = "polit_sesrq", data = d1, weights = NULL,
                     Ylabel = "Response Time (Logged)", Dlabel = "Presidential Co-Partisan", Xlabel="Politicization Ratio", na.rm=T)
# output
#pdf("OUTPUT FILE PATH",height=5,width=8)
#ff1
#dev.off()

# FIGURE G2: Binning Estimator Plot
ff2a <- inter.binning(Y = "rt_log", D = "presparty", X = "polit_sesrq", Z = vars, data = d1, nbins = 3,
                      Ylabel = "Response", Dlabel = "Presidential Co-Partisan", Xlabel="Politicization Ratio",na.rm=T)
ff2b <- inter.binning(Y = "rt_log", D = "presparty", X = "polit_sesrq", Z = vars, data = d1, nbins = 5,
                      Ylabel = "Response", Dlabel = "Presidential Co-Partisan", Xlabel="Politicization Ratio",na.rm=T)


    